function cheb_test(A,e,v)
% chebyschev acceleration of an eigenvalue problem

%clear;
tmp = 0;
if (tmp==0) % spd, e's > 0
    N = 10;
    A = diag( 2*ones(N,1), 0 ) + ...
        diag( -ones(N-1,1), 1 ) + ...
        diag( -ones(N-1,1),-1 ) ;
    [v,e] = eig(A); e = diag(e); e = sort(e, 'descend')
    v_ref = v(:,1);
elseif (tmp==1) % e's real
    N = 100;
    A = diag( ones(N,1)/6, 0 ) + ...
    diag( -5*ones(N-1,1),1) + ...
    diag( -5*ones(N-1,1),-1) ;  A(1,1) = 1; A(1,2) = -1; 
    A = A/10.4;
    [v,e] = eig(A); e = diag(e);
    v_ref = v(:,end);
elseif (tmp==2) % e's complex, but small im part (like MR)
    N = 6;
    A = [ ...
   -4.2375   -0.4520    0.6688    0.7525   -0.2964   -0.4013
    0.2553   -4.1710    0.4854    0.8391   -0.6815   -0.8261
    0.4117   -0.9044   -2.8905    0.8439    0.9657   -0.8898
    0.8897    0.2785   -0.5265    3.8175   -0.6881    0.4641
   -0.4014    0.8175   -0.5276    0.0077    3.4562    0.6913
   -0.7859    0.0262    0.9686    0.0141    0.4218    4.2588 ];
    %A = A/10.4;
    [v,e] = eig(A); e = diag(e);
    v_ref = v(:,end);  
elseif (tmp==3)
    N = 12;
    A = [ ...
   -4.2375   -0.4520    0.6688    0.7525   -0.2964   -0.4013
    0.2553   -4.1710    0.4854    0.8391   -0.6815   -0.8261
    0.4117   -0.9044   -2.8905    0.8439    0.9657   -0.8898
    0.8897    0.2785   -0.5265    3.8175   -0.6881    0.4641
   -0.4014    0.8175   -0.5276    0.0077    3.4562    0.6913
   -0.7859    0.0262    0.9686    0.0141    0.4218    4.2588 ];
    %A = A/10.4;
    A = [A A; A A];
    [v,e] = eig(A); e = diag(e);
    v_ref = v(:,end);  
else
    load MRb; 
    A       = MR;
    N       = length(A);
    [v,e]   = eigs(A,2);  % interesting, eigs(A,2) picks up 1&3, *not* 1&2
    e       = diag(e);  
    v_ref   = v(:,1);
    % need to load dat_A_v_e
    %N=length(e);
end
lambda_ref = e(1);

if (lambda_ref < max(e))
    disp(' warning: e(1) not max ev ' )
end
boo        = sort(abs(e));
rho_ref    = boo(end-1)/boo(end);
dom_idx    = find(e==max(boo));  % index of dominant m
v_ref      = v(:,dom_idx)/norm(v(:,dom_idx),2);
lambda_ref = e(dom_idx);
m     = sprintf('reference ev = %2.10f ', lambda_ref );
disp(m)


Vstore = 2*ones(N,1)+1; 
Vstore  = Vstore/norm(Vstore,1);
lo      = 1;  
errl    = 1; 
l_1     = 1; 
errV    = 1; 
l_2     = 0;
it      = 0;
p       = 1;
V_1     = Vstore; 
V_0     = V_1;
it2 = 0; rhoB = 1;

% save lala2
% return
%    0.999999922699637
%    0.994034661598560
%   -0.992856415680385
% load lala2
% eigs(A,3)          pi: 1329 
% return
% tic
tic
mxit = 1000;
nPi = 10; nCi = 100; rho = rho_ref;
while ( (errl > 1e-8 || errv > 1e-6 ) && it < mxit )
%     if (it<=100)
%         nPi = 30;
%     else
%         nPi = 10;
%     end
    if it2 < nPi % do standard power iteration
        V_2 = A*V_1/l_1;
        l_2 = l_1*norm(V_2,1)/norm(V_1,1);
        it2 = it2+1;  
    else
        [a,b] = alpha(p,rho);
        V_2   = A*V_1/l_1;
        V_2   = V_1 + a * ( V_2 - V_1 ) + b * ( V_1 - V_0 ); 
        l_2   = l_1*norm(A*V_2,1)/norm(A*V_1,1);
        p     = p + 1;
    end
    if (p==nCi)
        p = 1;
        it2 = 0;
    end
    errVo = errV;
    errV  = max( abs( V_2-V_1 ) ); % might not want abs!!!
   % if (it2==nPi && p==1)
   %     rho = rho_ref;
   %     rho = errV/errVo;
   %    % disp([' fuck = ', num2str(rho) ])    
   % end
    
    R2    = V_2-V_1; 
    R1    = V_1-V_0;
    V_0   = V_1;
    V_1   = V_2;  
    l_1   = l_2;    
    %rho   = 0.939658145053100;
    if (it2==nPi)
        rhoB  = sqrt( (R2'*R2) / (R1'*R1) );
    end
    % disp([' fuckB = ', num2str(rhoB) ])
    errl  = abs(l_2-lambda_ref);
    errv  = max( abs( V_2/norm(V_2)*sign(V_2(1)) - v_ref*sign(v_ref(1)) ) );
    
    
    it    = it + 1;
    if ( mod(it,10)==0 )
        m     = sprintf('      ev = %2.10f (%2.10f, %2.10f) rho = %2.6f in %i its',...
            l_2,errl,errv,rho,it );
        disp(m)
    end
end
m     = sprintf('final ev = %2.10f (%2.10f, %2.10f) rho = %2.6f in %i its',...
    l_2,errl,errv,rho,it );
disp(m)
toc
     
end

function [a,b] = alpha(p,rho)
        bb = 0;% -0.992856415680385/0.999999922699637; % we expect the most negative eigenvalue to have a 
                   % magnitude similar to the second eigenvalue (if it's
                   % not itself the second eigenvalue w/r to modulus).
    if ( p == 1 )
       % a = 2/(2-rho);
        a = 2/(2-rho-bb);
        b = 0;
    else
%         y = 2/rho - 1;
        y = (2-rho-bb)/(rho-bb);
%         a = 4/rho * ( Tn(p-1,y)/Tn(p,y) );
        a = 4/(rho-bb) * ( Tn(p-1,y)/Tn(p,y) );
%         b = (1-rho/2)*a - 1;
        b = Tn(p-2,y)/Tn(p,y);
    end
end


function z=Tn(p,y)
    z=cosh(p*acosh(y));
end
